/* Copyright 2012 Tobias Marschall
 *
 * This file is part of CLEVER.
 *
 * CLEVER is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CLEVER is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CLEVER.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <bamtools/api/BamAux.h>

#include "MismatchWeightTrack.h"

using namespace std;

MismatchWeightTrack::MismatchWeightTrack() {}

MismatchWeightTrack::~MismatchWeightTrack() {
	all_tracks_t::const_iterator it = tracks.begin();
	for (; it != tracks.end(); ++it) {
		delete it->second;
	}
}

void MismatchWeightTrack::extractFromAlignment(const BamTools::BamAlignment& aln, double alignment_weight, int phred_offset, std::vector<MismatchWeightTrack::mismatch_weight_t>* target) {
	assert(target != 0);
	if (!aln.IsMapped()) return;
	unsigned int ref_pos = aln.Position;
	unsigned int query_pos = 0;
	for (size_t i=0; i<aln.CigarData.size(); ++i) {
		const BamTools::CigarOp& c = aln.CigarData[i];
		switch (c.Type) {
			case 'M':
				assert(false);
			case 'S':
			case 'I':
				query_pos += c.Length;
				break;
			case 'N':
			case 'D':
				ref_pos += c.Length;
				break;
			case '=':
				ref_pos += c.Length;
				query_pos += c.Length;
				break;
			case 'X':
				for (size_t j=0; j<c.Length; ++j) {
					double w = alignment_weight * (1.0 - pow(10.0,-((double)(aln.Qualities[query_pos] - phred_offset))/10.0));
					assert(w >= 0.0);
					int w_int = (int)round(w * 10.0);
					if (w_int > numeric_limits<unsigned char>::max()) {
						w_int = numeric_limits<unsigned char>::max();
					}
					// cerr << "Adding " << w_int << " to position " << ref_pos << endl;
					target->push_back(mismatch_weight_t(aln.RefID, ref_pos, (unsigned char)w_int));
					ref_pos += 1;
					query_pos += 1;
				}
				break;
			default:
				break;
		}
	}
}

void MismatchWeightTrack::addAll(const std::vector<mismatch_weight_t>& mismatch_weights) {
	track_t* track = 0;
	for (size_t i=0; i<mismatch_weights.size(); ++i) {
		const mismatch_weight_t& m = mismatch_weights[i];
		if ((i==0) || (m.chromosome_id != mismatch_weights[i-1].chromosome_id)) {
			all_tracks_t::iterator it = tracks.find(m.chromosome_id);
			if (it == tracks.end()) {
				track = new track_t(m.position + 1000, 0);
				tracks[m.chromosome_id] = track;
			} else {
				track = it->second;
			}
		}
		if (m.position >= track->size()) {
			track->resize(m.position + 1000, 0);
		}
		assert(m.position < track->size());
		int w = track->at(m.position);
		w += m.weight;
		if (w > numeric_limits<unsigned char>::max()) {
			w = numeric_limits<unsigned char>::max();
		}
		track->at(m.position) = (unsigned char)w;
	}
}

std::auto_ptr<std::vector<MismatchWeightTrack::snp_t> > MismatchWeightTrack::getSnpCandidates(double weight_cutoff) {
	auto_ptr<vector<snp_t> > result(new vector<snp_t>());
	all_tracks_t::const_iterator it = tracks.begin();
	for (; it != tracks.end(); ++it) {
		for (size_t pos=0; pos<it->second->size(); ++pos) {
			double w = ((double)(it->second->at(pos))) / 10.0;
			if (w >= weight_cutoff) {
				result->push_back(snp_t(it->first, pos, w));
			}
		}
	}
	return result;
}
